function [ Fgrid ] = Matrix_Bivariate_4pointRefinement_boundary( grid , omega )
%-----------------------------------------------------------------
% Input:
% grid(m,n) - a grid over Z^2
% Output:
% Fgrid(2*m-1,2*n-1)
%
% ASSUMPTION: m,n > 4
%-----------------------------------------------------------------
% ABSTRACT
% Tensor product, includes boundaries. the middle point is insert in the order y->x of
% calculation
%-----------------------------------------------------------------
% NIr Sharon, 26-05-12
%-----------------------------------------------------------------

m = size(grid,1);
n = size(grid,2);
N = size(grid,3);

Fgrid = zeros(2*m-5,2*n-5,N,N);


%-----------------------------------------------------------
% First main loops: For vertices (new and old) on the grid

for x=1:(m)
    for y=1:(n)
        % Interpolation
        Fgrid(2*x-1,2*y-1,:,:) = squeeze(grid(x,y,:,:));
        % Y direction
        if (y~=n)
        if (y==1)
            Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary(squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)) , squeeze(grid(x,y+3,:,:)),omega);
        else
            if (y==(n-1))                       
               Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage_boundary(squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y-1,:,:)) , squeeze(grid(x,y-2,:,:)),omega);       
            else           
               Fgrid(2*x-1,2*y,:,:) = Matrix_FourPointAvarage(squeeze(grid(x,y-1,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x,y+1,:,:)) , squeeze(grid(x,y+2,:,:)),omega);                    
            end
        end
        end
        % X direction
        if (x~=m)
        if (x==1)        
            Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary(squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)) , squeeze(grid(x+3,y,:,:)),omega);
        else
            if (x==(m-1))            
               Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage_boundary(squeeze(grid(x+1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x-1,y,:,:)) , squeeze(grid(x-2,y,:,:)),omega);
            else
               Fgrid(2*x,2*y-1,:,:) = Matrix_FourPointAvarage(squeeze(grid(x-1,y,:,:)) , squeeze(grid(x,y,:,:)) , squeeze(grid(x+1,y,:,:)) , squeeze(grid(x+2,y,:,:)),omega);
            end
        end
        end
    end
end

%-----------------------------------------------------------
% Second main loops : grid face centers 

for x=1:(m-1)
    for y=1:(n-1)
        % Y direction        
        if (y==1)      
            m_y = Matrix_FourPointAvarage_boundary(squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y+3,:,:)) , squeeze(Fgrid(2*x,2*y+5,:,:)),omega);
        else
            if (y==(n-1))            
                m_y = Matrix_FourPointAvarage_boundary(squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y-3,:,:)) , squeeze(Fgrid(2*x,2*y-5,:,:)),omega);
            else
                m_y = Matrix_FourPointAvarage(squeeze(Fgrid(2*x,2*y-3,:,:)) , squeeze(Fgrid(2*x,2*y-1,:,:)) , squeeze(Fgrid(2*x,2*y+1,:,:)) , squeeze(Fgrid(2*x,2*y+3,:,:)),omega);
            end
        end
        % X direction
        if (x==1)        
           m_x = Matrix_FourPointAvarage_boundary(squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x+3,2*y,:,:)) , squeeze(Fgrid(2*x+5,2*y,:,:)),omega);
        else
            if (x==(m-1))
                m_x = Matrix_FourPointAvarage_boundary(squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x-3,2*y,:,:)) , squeeze(Fgrid(2*x-5,2*y,:,:)),omega);
            else
                m_x = Matrix_FourPointAvarage(squeeze(Fgrid(2*x-3,2*y,:,:)) , squeeze(Fgrid(2*x-1,2*y,:,:)) , squeeze(Fgrid(2*x+1,2*y,:,:)) , squeeze(Fgrid(2*x+3,2*y,:,:)),omega);
            end
        end
       % dif(x,y) = m_x-m_y;
        Fgrid(2*x,2*y,:,:) = GM2(m_x,m_y,.5);
    end
end
  
end



